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In this paper we calculate and visualize the dynamics of an ensemble of electrons trapping in an electrostatic 
wave of slowly increasing amplitude, illustrating that, despite disordering of particles in angle during the trapping 
transition as they pass close to X-points, there is still an adiabatic invariant for the great majority of particles that 
allows the long-time distribution function to be predicted. Possible application of this approach to recent work 
on the nonlinear frequency shift of a driven wave is briefly discussed. 
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1 Introduction 

The Bernstein-Greene-Kruskal (BGK) [1] Vlasov con- 
struction shows that fully nonlinear collisionless electro- 
static waves and shocks, time-independent in some frame 
of reference, can in principle exist, but makes no statement 
as to the physical accessibility of these structures. The 
physical question is, given an initially time-independent, 
spatially uniform distribution function, which evolves into 
a nonuniform state either due to instability or external forc- 
ing, what is the final long-time steady state after coarse 
graining/phase-mixing the distribution function? (Also, in 
which frame is it time-independent?) 

Consider a charged particle moving in a nonlinear 
wave propagating in the z-direction with electrostatic po- 
tential in the laboratory frame given by 



<f>(z,t) = u{9\et) , 



(1) 



where u{6) is 2^-periodic, u(6 + 2n) = u{6). In many cases 
u can also be assumed to have half-period antisymmetry, 
u(0 + n) = -u{6) (i.e. it has only odd harmonics) and we 
shall assume this. The et dependence expresses the possi- 
bility that the amplitude and waveform can evolve slowly 
with time, the adiabatic limit being defined by e — > 0. 
The phase angle 9 is defined by 
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e = kz-\ (o(et')dt' , 



(2) 



where the wave vector ke z is constant, but we have allowed 
for a) to evolve with time, e.g. due to a nonlinear frequency 
shift. The wave phase 9 at the particle position forms a con- 
venient generalized coordinate for describing the particle 
dynamics as it is nondimensional. Also, 9/k = v z - v P h is 
the particle velocity in the wave frame, i.e. after a Galilean 
transformation to a frame moving at the phase velocity 
v p h = co/k. (Note that the wave frame is noninertial if oj is 



time-dependent.) The equation of motion for a particle of 
mass m, charge q is 
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(3) 
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The particle dynamics during the evolution between 
the initial and final states cannot be solved analytically be- 
cause the energy is not a constant of the motion. How- 
ever, Dewar [2, 3] showed that the approximation of adia- 
batic invariance, combined with phase mixing, provided a 
sufficiently accurate mapping between the initial and final 
distribution functions that reasonable estimates of nonlin- 
ear frequency shift and saturation amplitude of an unstable 
Langmuir wave could be made analytically. 

Motivated by applications in laser fusion and beam 
physics, the utility of the adiabatic approximation has re- 
cently been confirmed numerically by Lindberg et al. [4] 
for a driven Langmuir wave. In [4] the long-time wave 
response was calculated both by particle simulation and 
semi-analytically using the adiabatic Vlasov approxima- 
tion [2,5] or the electron dynamics in a wave field of slowly 
varying amplitude and frequency. The authors allowed for 
a nonsinusoidal waveform due to excitation of harmonics 
at large amplitudes, and for the generation of a dc electric 
field to satisfy the assumed external circuit conditions of 
zero initial and final spatial-mean current (i.e. the dc com- 
ponent of the current). Both methods of calculation agreed 
well, confirming the utility of the adiabatic approximation. 

In the small-amplitude limit, where the waveform is 
sinusoidal, Lindberg et al. [4] also calculated a nonlinear 
frequency shift analytically in terms of elliptic integrals, 
finding an expression proportional to the square root of the 
final wave amplitude, with a coefficient in agreement with 
the earlier calculation by Dewar [2]. However, Fig. 4 of [4] 
shows quite poor agreement between this small-amplitude 
formula and their finite-amplitude simulations and semi- 
analytic calculations, except at very low amplitudes. 
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In Sec. [2] we specialize to a sinusoidal wave of fixed 
frequency and describe a numerical experiment in which 
a set of electrons of the same initial wave-frame energy, 
but with different phases, is evolved as the amplitude is in- 
creased. A special choice of growth function allows the 
extreme adiabatic limit to be probed and it is shown that 
adiabatic theory is statistically very accurate in this limit, 
even though particles originally launched within one wave- 
length of each other end up being trapped in several wave 
troughs. This is investigated in more detail in Sec.[3]where 
the details of the trapping process are visualized, illus- 
trating that close encounters with hyperbolic X-points can 
break up the ordering of the particles. 

Section |4] shows the result of evolving an ensemble 
of particles, some of which are trapped while others re- 
main free, illustrating why the coarse-grained distribution 
function of trapped particles is half its initial pre-trapping 
value. This is used in Sec. [5] to explain the basis of our 
earlier formalism [2, 3]. 

In the concluding remarks in Sec. [6] we mention some 
other early theoretical work [6, 7] that may point the way 
to improving the agreement between the asymptotic am- 
plitude expansion for the nonlinear frequency shift and the 
simulation results of Lindberg et al. through the inclusion 
of higher order terms. 

2 Adiabatic trapping 

To understand adiabatic trapping more clearly, in this sec- 
tion we study a specific example, a set of electrons with 
charge q — -e moving in a sinusoidal wave 

0(z,f) = 0i(Ocos0. (4) 

We here assume co constant, so 9 = kz - cot. 
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Fig. 1 Adiabatically increasing amplitude of the electron poten- 
tial energy, ei/>, in the wave vs. time t, between t = -t f 
and t = tf as described in the text. 

In the autonomous case, e<p\ = const, the wave-frame 
energy 

W=^6 2 + q4> (5) 



is a constant of the motion, which allows solution of this 
one-degree-of-freedom system by quadratures. For the si- 
nusoidal wave there is a direct analogy with the physical 
pendulum, where 9 is the angle relative to the vertical. The 
motion of an autonomous nonlinear pendulum is known to 
be solvable in Jacobian elliptic functions [8]. For a pen- 
dulum of length Z, the analogue of ek 2 cb\jm is g/l. Thus, 
for the pendulum the analogue of a changing wave ampli- 
tude would be a changing gravitational field g(t). Trapping 
in a wave is analogous to the transition from rotation to 
libration for a pendulum. 

In this section we investigate the dynamics of a set of 
initially untrapped electrons moving in a wave potential of 
the form specified in Eq. Q, with the amplitude function 
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and with cf> 9 = 1.5 and Acf> — 1.2 in units such that e — k — 
m = 1. This function, depicted in Fig.[T] is chosen so that 
its time rate of increase is zero at the beginning, t = -tf, 
middle, t = 0, and end, t = tf, of the calculation period. 
This choice makes the adiabatic description very accurate 
while keeping the total time of the calculation reasonably 
short: defining tf = K/ecJ^, where = kiefi^/m) 1 ? 2 is the 
bounce frequency at the bottom of a wave potential trough 
at time t — 0, the results in this paper were obtained with 
e = 0.04, giving t f = 64.13. 
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Fig. 2 Phase-space positions of the set of 101 particles described 
in the text at the final time t = t f . 

Figure|2]shows the final positions of a set of 101 elec- 
trons initialized at t = -tf with the same initial total energy, 
W = 1.22515, but at different initial points (6, 0) on the up- 
per W-contour, i.e. with positive velocity 6. The particles 
were initialized in an interval of 9 of width 27r whose end- 
points were chosen to the left of the origin, such that all 
the particles trapped in wave troughs near the origin. Each 
electron was assigned a unique color, starting from red at 
the left initial endpoint to blue on the right initial endpoint. 

Rather than using action-angle variables [4] we adopt 
a more Lagrangian approach and study the particle orbits in 
9, 9 phase space. It is a standard result of adiabatic invari- 
ant theory that the phase-space area Jfd9d9 under the con- 
tour W = C(et), where C denotes a constant of the motion 
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in the autonomous case, is an adiabatic invariant for both 
trapped and passing (free) particles. However, Best [9 J] 
discovered that this is true for most particles even through 
the trapping or detrapping transition, provided the 6 in- 
tegral is taken over one 2n period between the upper W- 
contour and the 6>-axis while the particle is free, or over 
half a W-contour when the particle is trapped. 

Following Dewar [2, 3] we divide the phase-space 
area by 2n, so that the adiabatic invariant becomes the 8- 
average of the positive branch of the solution for the angu- 
lar velocity of the equation W = const, 



0(W) = 
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)(W-q<f>) l/2 (7) 



where H(-) is the Heaviside step function. 

Equation (|7]) is valid as an adiabatic invariant for an 
arbitrary waveform as in Eq. ([T} provided the nonlinear 
frequency shift can be ignored. In the special case of a 
sinusoidal wave, the integration can be done analytically 
in terms of complete elliptic integrals E(m) and K(m) [8]. 
For passing particles we have 
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while for trapped particles 
4w b 



G(W) 



[E(w)-(l -w)K(w)] 



(8) 



(9) 



where w(t,W) = (W + <?0i)/(2e^i) is > 1 for passing 
particles and < 1 for trapped particles. Here Wb(0 = 
k(e<pilm) l l 2 is the bounce frequency at the bottom of a 
wave trough at time t . 

In Fig. [2] the final separatrix, i.e. the contour W = 
e(f>\(tf) passing through the X-points at odd multiples of 
7r, is shown in yellow. Trapping comes about due to the 
expansion of the separatrix as e(p\ increases — when the 
area enclosed by the separatrix exceeds the area under one 
period of the upper initial W-contour the particles must 
trap to conserve the adiabatic invariant. This was used to 
choose the initial value of W so that the adiabatic theory 
prediction for the trapping time was t = 0. 

The adiabatic prediction for the final energy is shown 
in green, and it is seen that, at the resolution of the plot, all 
the particles do lie on the contour predicted by Best's the- 
ory, but the higher resolution W-0 representation in Fig. [3] 
shows departures from adiabaticity. However, even the 
least adiabatic points are quite close to adiabatic and most 
points cluster close to the adiabatic prediction. 

In Fig. [4] we quantify this observation of clustering, 
verifying that the great majority of particles do lie close 
to the adiabatic prediction. The mean energy at t = +tf 
is W = 1.7072 with a standard deviation 0.002, which is 
0.1% of the mean. However, the standard deviation greatly 
overestimates the discrepancy between W and the adiabatic 
prediction, W ad = 1.70687, which is within 0.02% of the 
mean. 
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Fig. 3 Particles from the initially well-ordered set are, at / = 
+tf, distributed over three different wave troughs but have 
total energies very close to the adiabatic prediction. 
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Fig. 4 Histogram showing the statistical distribution of final en- 
ergies of trapped particles. 
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Fig. 5 



Onset of trapping in angle space for a subset of 1 1 parti- 
cles starting on an energy contour outside the separatrix 
with a spread of 2rc in angle, but which trap in three dif- 
ferent wave troughs. 



See also Elskens and Escande [10] and references therein. 
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3 Details of trapping 

We see in Figs. [2] and [3] that the initial ordering in 8 has 
been severely disrupted by the final time, with the parti- 
cles, initially within the same 2n interval in 8, now spread 
over three different wave troughs. To visualize the mech- 
anism for this disordering we now examine the trapping 
transitions more carefully. 




Fig. 6 Angular velocity for two adjacent orbits from the set used 
in Fig. [5] which trap in different wave troughs. The mean 
angular velocity evolves slowly until rapidly dropping to 
zero at the onset of trapping. 



As is seen in Figs.[5]and[6] the particles do indeed pass 
through the separatrix and become trapped in the vicinity 
of t — 0, when the rate of increase of the amplitude is very 
small. Because of this, the effective e is much smaller dur- 
ing the trapping period (when adiabatic theory is least ac- 
curate) than its actual value of 0.04. Comparing the slope 
calculated at a typical trapping time (taken to be ~ 5 on 
the basis of the results shown in Fig. [5]), of a simpler ramp- 
up function, <f>® + A(f>smnt/2tf, we estimate the effective 
adiabatic expansion parameter to be e e ff ~ 2 x 10" 3 . 

It is also instructive to examine the trapping process in 
phase space. We see in Figs.[7]and[8]that all particles come 
very close to the hyperbolically unstable X-points before 
they finally trap, thus amplifying the small differences in 
the individual histories of each orbit and accounting qual- 
itatively for the disordering in 8 seen at the final time in 
Figs. [2] and [3] We now see that this disordering occurs 
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7 Phase-space plots of the trapping of the set used in Fig. [5] 
showing the two trapping troughs to the right of the 6 
origin. 
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Fig. 9 Initial state of ensemble, at t = —tf. 

during trapping and is presumably not primarily due to the 
accumulation of phase differences during the subsequent 
spiraling in toward the final adiabatic energy, though the 
small deviations from this energy seen in Fig. [8] will con- 
tribute to disordering from small variations in the nonlinear 
bounce frequencies. 

4 Trapping of an ensemble 

In Fig. [9] we show an ensemble of initial points (in gray) 
randomly distributed with equal phase-space probability 
density, f(6, 8, -tf) = const, between two W-contours. The 
lower contour is the same one used in Sec. [2] above, with 
the same color-coded set of 101 initial points shown, but 
projected modulo 2tt onto the interval -n, n. That is, we 
are now using the cylindrical topology of the pendulum 
phase space (or, alternatively, duplicating the initial set of 




Fig. 10 Final state of ensemble, at t = +tf. 
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points across all periods of the wave). The upper contour 
is chosen so that electrons on it never trap, and 101 color- 
coded initial points on this contour are also shown. 

Figure [10] shows the final state after the points shown 
in Fig. [9] are evolved dynamically from t = -tf to t — tf. 
The trapped points on the inner adiabatic contour (green) 
are the same ones shown in Fig. |2]but now projected into 
one period as explained above. The un trapped points out- 
side the separatrix on the upper contour have evolved adi- 
abatically to a very high accuracy as they never come near 
the expanding separatrix (yellow). 

As the area of phase-space elements is preserved by 
the dynamics, the probability density in the untrapped re- 
gion, between the upper separatrix and the upper contour, 
is the same as used initially, / pass (0, 6, tf) = f(9, 6, -tf). 
However the trapped points, between the inner contour and 
the entire separatrix, are now distributed quasi-randomly 
over twice the area they occupied initially (because the 
adiabatically conserved area is that of only the inner half- 
contour above the 0-axis). Thus the coarse-grained prob- 
ability density of the trapped particles is half the initial 
value: 

fu^{0,e,tf) = \f{e,e,-tf). (io) 

5 Distribution function after adiabatic 
wave excitation 

Here we relate the above result to the formalism used in 
our earlier calculations [2, 3]. We work in z, v z phase space 
in the lab frame, so the relation to the wave frame must be 
made explicit. Thus Eq. Q is now written 

1 7 

W = -m(v z - v ph ) z + qcf>(z - v ph r, st) , (11) 

where v z = z and phase velocity v p h may include a nonlin- 
ear frequency shift. Solving Eq. (|TTJ for v z we have two 
solutions, corresponding to particles going to the right or 
left in the wave frame, v* = v p h ± m, where 

u = (2/m) l/2 [W - q<p(z - v ph t,st)] l/2 . (12) 

We assume there to be no wave at t = -oo, so then u 
reduces to uq(W) = (2/m)^ 2 VK^ 2 and the initial distribu- 
tion function /o(v z ) can be written as a two-branched func- 
tion of W, which we denote by 

Ft oa (W) = Mv ph ±u (W)). (13) 

The long-time, coarse-grained (or phase averaged [11]) 
distribution function / ra (z, v z , t) will be independent of time 
in the wave frame, and so must also be a function of the 
constant of the motion W, which we denote by F X (W): 

ft»(W)s/fc±|vJ,f = +«) (14) 

The goal of adiabatic theory is to find the transformation 
between F_ M (W0 (or f ) and F +0O (W). 



Restricting to the case where the nonlinear frequency 
shift and dc field can be ignored, we use the adiabatic in- 
variant 6(W) defined in Eq. {7}, or, rather, u(W) = 0(W)/k. 
We also see from the previous section that a transition from 
a passing orbit at t = -oo to a passing orbit at t — +oo pre- 
serves both the value of the distribution function and the 
direction of the wave-frame velocity (provided there is no 
intermediate trapping), while a transition from a passing 
orbit to a trapped orbit mixes initially left-going and right- 
going passing particles in the same band of W and halves 
the value of their individual distribution functions, 

Thus the transformation between /q and F +00 is as 
summarized below. 

• Passing to passing: 

FtJW) = Mv ph ±u(W)) (15) 

• Passing to trapped: 

F^iW) = i[/o(Vph-«(W)) 

+ /o(v„h + «(W))] . (16) 

6 Conclusion 

The numerical and graphical study presented here verifies 
the applicability of the adiabatic approximation to the cal- 
culation of the long-time coarse-grained distribution func- 
tion of a plasma after the growth and saturation of a slowly 
growing instability or driven wave. Yet to be studied in 
similar detail is the effect of an external electric field, and 
of a frequency that changes with time, in order to verify 
the adiabatic theories presented in previous work [3,4]. 

Also awaiting further work is a careful matched 
asymptotic expansion using adiabatic theory to calculate 
the nonlinear frequency shift of a driven wave beyond the 

1 /2 • 

first Oitpy ) term in order to compare with the numeri- 

1 /2 

cal calculations of Lindberg et al. [4]. While the 0((f>^ ) 
term due to trapped particles must dominate for very small 
amplitudes, higher powers of (f>\ will become dominant at 
higher amplitudes, and it is clear this must be occurring for 
the parameters used in Fig. 4 of [4]. 

The 0{<p 2 ) frequency shift due to the nonlinear re- 
sponse of the bulk of the distribution function, which is not 
subject to particle trapping effects, was calculated earlier 
using a "waterbag" distribution function (effectively a fluid 
model) by Dewar and Lindl [6] and using an averaged- 
Lagrangian oscillation-center kinetic method by Dewar 
[7]. Winjum et al. [12] have recently postulated that such 
a "fluid" nonlinear frequency shift can simply be super- 

1 11 

imposed on the 0{<p^ ) frequency shift from trapped par- 
ticles, but it seems a priori quite possible that a careful 
asymptotic expansion will reveal terms at intermediate or- 
ders, <9(0i ) and 0((p\ /2 ). 
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